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Abstract. - We uncover a solvable generalization of the Kuramoto model in which shears (or 
nonisochronicities) and natural frequencies are distributed and statistically dependent. We show 
that the strength and sign of this dependence greatly alter synchronization and yield qualitatively 
different phase diagrams. The Ott-Antonsen ansatz allows us to obtain analytical results for 
a specific family of joint distributions. We also derive, using linear stability analysis, general 
formulae for the stability border of incoherence. 



Introduction. — Collective synchronization is a com- 
monly observed phenomenon in nature and in some tech- 
nological applications [IrE], in which mutual interactions 
succeed to entrain the rhythms of a heterogeneous ensem- 
ble of self-sustained oscillators. It is mathematically cap- 
tured by a prototypic minimal model put forward by Ku- 
ramoto more than thirty years ago [2][6]- Moreover, this 
model is a suitable framework for the quantitative analysis 
1 of a variety of physical systems such as arrays of Josephson 
junctions [7J or mechanical rotors/oscillators [8TU0j. 

The universal form of a limit-cycle close to a Hopf bi- 
furcation led Kuramoto to analyse collective synchroniza- 
tion resorting to the mean-field version of the complex 
Ginzburg-Landau equation with disorder [5J[Tni[H] : 

K N 

Zj =z J [l + i{uJ ] +q )-{l + iq ] )\z J \ 2 ] + —^2{z k ~z 1 ), (1) 

fc=i 

where Zj = Qje j , and j — 1, . . . , N 3> 1. Here, uij is the 
natural frequency of the j-th oscillator, whereas qj is the 
so-called shear (or nonisochronicity) that quantifies the 
dependence of the oscillation frequency on the amplitude. 
Under the assumptions that the coupling is purely diffu- 
sive (K real) and weak (\K\ small), a phase reduction of 
eq. (|T]| yields 

K N 

9j = Uj + K qj + jjYl I sin ^ fc - 9 i) ~ q J C08 ( 9k ~ • ( 2 ) 

k=l 

We may also cast eq. ([2]) in a more compact form: 



with tan/3j = qj and \f3j\ < -|. Under the simplifying 
assumption that the shears are not distributed, qj = q 
(=> j3j — $), the so-called Sakaguchi-Kuramoto model [12] 



is recovered (redefining cj' = 



K tan $, and K' 



K tan Pj 



cos Pj N ^ 



1^-9,. + ft) (3) 



k=l 



Kj cos/3). Additionally, under the more severe constraint 
qj = 0, eq. becomes the standard Kuramoto model [B]. 

The goal of this work is to perform a detailed analysis 
of phase equations ^ under the assumption that the nat- 
ural frequency and the shear of each oscillator are drawn 
from a joint probability density function (PDF), p(u>,q). 
A particular case of this problem has been recently anal- 
ysed assuming the parameters uj and q to be independent 
random variables, p{ui,q) — g(ui)h(q) |13j . An interesting 
finding is that, if the width of the distribution h(q) ex- 
ceeds a precise threshold, diffusive coupling is unable to 
counteract shear heterogeneity leading to a self-organized, 
synchronous state. This result is in sharp contrast with 
the well-known prediction of the Sakaguchi-Kuramoto — 
or the Kuramoto — model, where collective synchroniza- 
tion is assured at large enough K values. 

How do these results translate into the case where nat- 
ural frequencies and shears are statistically dependent? 
This is the case one may encounter when studying the 
synchronization of any particular class of self-sustained os- 
cillators. Generally, model-specific parameters affect both 
the oscillator's natural frequency and shear. Therefore, 
heterogeneity in certain parameters will also result into 
heterogeneity of ui and q with some functional or statisti- 
cal dependence between them; see e.g., eqs. (7) and (8) in 
[7] for such a situation, though the heterogeneity of q is 
eventually neglected to simplify the analysis. In previous 
work, parameter dependencies were found to influence the 
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effect of diffusive coupling on the variance of the ensem- 
ble's oscillator frequencies, a phenomenon called 'anoma- 
lous phase synchronization' |14) . 

In this Letter, we define a conditional probability of uj 
given q, g c (ui\q), such that p(uj,q) = h(q)g c (uj\q). The 
marginal PDF h(q) is assumed to be unimodal, sym- 
metric and centred at qo. Additionally, the conditional 
probability g c is chosen to be unimodal and of the form 
g c (uj\q) = g c (u!—mq). This restricts our results to a partic- 
ular class of distributions that is nonetheless wide enough 
to illustrate a number of different synchronization scenar- 
ios (particularly depending on the sign of m). 

Continuum limit. — In our theoretical analysis we 
neglect finite-size effects and consider the thermodynamic 
limit N — > oo of model (|2|) . It is possible then to drop the 
indices and define the probability density for the phases 
/(#, uj, q, t). Thus f(9, uj, q, t) d9 duj dq is the ratio of oscil- 
lators at time t with phases between 9 and 9 + d9, natural 
frequencies between uj and uj + duj, and shear between q 
and q + dq. The density function / obeys the continuity 
equation 

dtf = -d e (L + Kq+§ h" i9 (! " " c - c -] } /) . 

(4) 

where c.c. stands for complex conjugate of the preceding 
term, and the complex order parameter r is 



r(t) = Re m = 



e ie f(0,u),q,t)d0du)dq. (5) 



The mean field r measures the degree of synchronization 
of the system. If the oscillators are uniformly distributed, 
a state commonly referred to as incoherence, f(9,uj,q,t) 
equals p(w, g)(27r)~ 1 , and r vanishes. States for which part 
or all of the population is entrained at a given frequency 
result in a nonuniform distribution of the phases such that 
R> 0. 

The density function f(9, uj, q, t) is real and 2-7r-periodic 
function in the 9 variable with the Fourier expansion 



2tt 



(6) 



where /; = /*;, fo = 1. Inserting this Fourier series 
into the continuity equation an infinite set of integro- 
diffcrcntial equations for the Fourier modes is obtained: 

Kl 

d t h = -U( u + Kq)fi + —[r*(l + iq)fi- 1 - r(l - iq)f l+1 ] 

(7) 

Note that the order parameter §5§ is only determined by 
the first Fourier mode: 



Ott-Antonsen ansatz. — Recently Ott and Anton- 
sen (OA) found an ansatz [TJ] , which is generically [TOUT?] 
satisfied by the asymptotic dynamics of the Kuramoto 
model (q = 0) — and, remarkably, of many variations of 
it, see e.g. [T51[T5H2Tj . In our case this ansatz takes the 
form 

f l (u,q,t)=a(u,q,t) 1 (9) 

This defines a family of solutions of eq. with the con- 
straint that a satisfies 

d t a = -i(LJ + Kq)a+—[r*(l+iq)-r(l- iq)a 2 } . (10) 

Our simulations indicate that the asymptotic solutions of 
the system indeed belong to the OA manifold. 

We consider first a family of joint PDFs p(u>,q) = 
h(q)g c (u)\q), with Lorentzian (Cauchy) marginal distribu- 
tion h: 

Hi) = 7 — ^ — o, (if) 

{q-qoY + i 2 

and Lorentzian conditional distribution g c : 



9c{u]\q) 



6/ir 



-u -m(q-qo)] +8 2 



(12) 



The specific family of PDFs defined by eqs. (fTTj) and (p~2|) 
allows us to obtain simple low-dimensional ordinary dif- 
ferential equations for the order parameter dynamics, and 
to tune the statistical dependence between w and q with 
the parameter m. (The case m = — u> and q indepen- 
dent random variables — was already addressed in [13].) 
In the limiting case 5 — > 0, g c becomes a Dirac's delta, 
and this results in a (deterministic) linear relationship: 
uj = ujq + m(q — qg). The terms r and r* in eq. (flTJ)) can be 
evaluated by means of the residue's theorem inserting the 
PDFs (fTTj) and (fT2|) in eq. |[8}, and closing the integration 
paths in the complex plane. Concerning variable uj, the in- 
tegration must be done in the lower half complex w-plane 
because a can be analytically continued in that region, 
as occurs in the Kuramoto model, see [15] for details. In 
partial fractions g c (ui\q) = (27ri) _1 {[a; — ujq — m(q — qo) — 
i<5] -1 — [uj — uiq — m(q — qo) + i(5] -1 }, and the integration 
over uj in eq. ([8]) involves only the value of a at the pole 
uj p = uj — m(q — q ) — iS, see [T3]. Hence, eq. $8$ becomes 
in this particular instance 



h(q) a(uj = uj p , q, t) dq. 



(13) 



To evaluate this integral over q we must proceed more 
cautiously to warrant that a can be analytically extended 
into the suitable half g-plane (q = q r + ift, with either 
qi > or qi < 0). Equation (fTO]) for a at (uj p , q) is: 

dtOL — — \\jJo — mqo + q(m + K) — iS\a 

+ y [r*(l + iq)-r(l-iq)a 2 ] (14) 



r*(t) 



p(uj, q)fi(uj, q, t) duj dq. (8) If a = |a|e ^ is analytic, it satisfies the Cauchy-Riemann 

> conditions, and this can be demonstrated to imply d qr \a\ + 
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d qi \a\ > 0. In consequence, the maximum of \a\ is nec- 
essarily located on the boundary (namely, on the integra- 
tion contour). Under the assumption that a is analytic 
at t = 0, analyticity will hold for all t > if a remains 
finite because a is the solution of the ordinary differential 
equation (fTDJ (see Theorem 8.4 in Chapt. 1 of 22 ). More- 
over we require \a\ < 1 everywhere in the selected half 
complex q-plane; otherwise the Fourier modes diverge, see 
eq. ((9|). After some algebra, one finds that on the real 
g-axis, eq. (fl4|) yields: 



d t \a\ = -S\a\ + y Re [rV'^l 



1 



(15) 



which gives dt\a\ = —5 < at |a| = 1. This implies that 
if \a\ < 1 at t = 0, this will hold for all t < 0. On the 
contour closing at infinity q = l^e" 5 with |g| — > oo, the 
dominant contributions (of order \q\) at \a\ — 1 give: 



dt\a\ = 



- K(l — i?cos%)] |g| sin$, 



(16) 



where \ = ip(q,t) — $>(t). If dt\a\ < is fulfilled in cither 
•d £ (0,7r) or ■& 6 (— 7r,0), we can safely choose that path 
for the contour closing in the integration of (TT3l) . The 
problem now is that if m ^ there are values of K in 
eq. (fT6|) where the desired relation <9t|cv| < cannot be 
fulfilled due to the "uncontrolled" angle \- Instead of ig- 
noring those parameter values we shall make the assump- 
tion that solutions in the range R < R x where <9t|a| < at 
| a | = 1 is fulfilled (choosing the appropriate half-plane), 
can be correctly studied within this framework. Thus, 
eq. (|16l) dictates that the analysable range of R is bounded 
by 

R x =mm(l,\l + m/K\), (17) 

what in particular implies that, in principle, the stability 
of incoherence (R = 0) can be always determined, save 
at K = -m (R x = 0). We must take e (— tt,0) if 
m + K > 0, and d G (0, 7r) if m + K < 0, for the closing 
of the integration contour in eq. (|13[) . Thus, the order 
parameter is determined by the value of a at the poles 



r*(t)=a(LO = LoP,q = qP,t), 



(18) 



with q p = qo — yy for m + K > 0, and q p = qa + for 
rn + K < 0. Equation (|18l) yields the relations R(t) = 
\a(uj p , q p , t)\ and = ijj(u p ,q p ,t), and hence it suffices 
to study eq. (TTU]) at (u>,q) = (uj p ,q p ). 

Recalling that q p = qa^fij, andw p = uj Q +m(q p —q n ) — iS, 
we obtain that the modulus and the phase of the order 
parameter (inside the OA manifold) obey Stuart-Landau 
equations: 



R 



-(fTm 7 +|(lT7)(l-fl 2 ) 
K 



* : CJ + yq (l-i? 2 ) 



R (19) 
(20) 



Remarkably, the radial dynamics does not depend on go, 
something that stems from the peculiarities (pointed out 
in [T3]) of the Lorentzian distribution. 




Fig. 1: Phase diagrams for positive dependence between cj and 
q: eq. (|12|) with m > 0. (a) Purely linear dependence 5 = 0. (b) 
S = 0.2m. The solid lines correspond to the loci of bifurcations 
given by the analytic formulas (|21l) . (|23[) . and K = rn (see text). 
The dotted line is a bound of the region of bistability given by 
eq. (|24|l . Dashed lines are obtained from numerical simulations 
with N = 2000 (a), 40000 (b) oscillators. Our numerics showed 
good agreement with boundaries H21|) and Q23p . data not shown. 
In the simulations we took {ujj, qj}j=i,...,N deterministically 
to represent p(cu,q) given by eqs. and (|12[) ; the selected 
parameters were m = 1 and luq = go = s. 



In the incoming paragraphs we present separately the 
cases of positive and negative m, as these two cases yield 
qualitatively different results. 



> 0). 



In this case, 



1 for K/m > — h, and R x < 1 for 



Positive dependence (to 
eq. (fT7|) implies i?> 

if /to < — h ■ In the latter region we cannot solve the prob- 
lem completely within the OA framework because possible 
attractors with R g [R x , 1] are not captured by the theory. 

UK > —to, the signs "=p" in eq. (|T9)) must be replaced 
by "— ". It can be easily seen that incoherence is stable 
everywhere, except above the line: 



K?> = 



2(to7 + 5) 
1-7 



with 7 < 1 



(21) 



A phase diagram for two values of 8 can be seen in fig. [TJ 
At a supercritical bifurcation gives rise to a partially 
synchronized solution with 



R 2 



K - Kg 

K 



(22) 



Remarkably, this formula coincides with the one obtained 
in the standard Kuramoto model [2] (recovered at 7 = 
qo = 0). 

If K < -to, one must replace "=f" by "+" in eq. (fl9|) . 
The resulting equation predicts the incoherence to be un- 
stable in the wedge-shaped region between K = —to and 



(2) _ 2(to 7 -<5) 
1+7 



with 7 > 1 



2rf 



(23) 



If ui and q are let to be progressively less statistically de- 
pendent (m — » 0), the tip of this region goes to 7 = 00. 
Thus, the interval of 7 where incoherence is stable for all 
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K becomes infinite as m — » 0, in consistence with our 
result in (T3] for the independent case (to = 0). Inside 
the wedge-like region where incoherence is unstable we 
can presume — and confirm numerically — the existence of 
a stable partially synchronized solution (with R > i? x ). 

(2) 

Moreover, as the instability of incoherence at K c is sub- 
critical, we can infer the existence of a region of bistability 
incoherence-synchronization below this line. The unstable 

,(2) 



solution with R > appearing at Kc 
determined up to 



can be analytically 



K h = - 



TO 2 (1 + 7 ) 

2(m + S) 



(24) 



where it acquires an R larger than R x . Hence Kb is a 
bound (surprisingly tight) for the region of bistability, see 
fig. 1. 

Negative dependence (m < 0). — In the case of 
negative to, eq. (fTTj) tells us that our eqs. (fTO|) and (|2"U|) 
apply to all R values when K/\m\ < I, while otherwise 
their validity only holds in a certain range R < R x . In 
contrast to the case of positive to, the phase diagram un- 
dergoes several transformations as the ratio between 6 and 
|m| varies. Next we describe the three main situations sep- 
arately, see fig. 2. 

Case I: < 5 < \m\/2; fig. 2(a,b). If K > \m\ inco- 
herence is stable only above the line 



2(-\m\~f + 8) 



1 



7 



with 7 > 1 



(25) 



where an unstable solution branches off incoherence obey- 
ing relation (f22j). As presumable, a region of bistability 
between incoherence and synchronization (with R > R x ) 
is found. For K < |m| incoherence is stable everywhere 
except above the line 



2(|to| 7 + <S) 
1+7 : 



with 7 < 1 - Mr 



(26) 



where it undergoes a supercritical bifurcation. 

Our numerical simulations reveal that a stable coherent 
solution exists below K/\m\ = 1 in the region of stable in- 
coherence, see bottom panels of fig.[2](a,b). For 6 — 0, this 
solution is continuation of a fully synchronized solution 
existing at K/\m\ = 1 with R = f™^ h(q)(l + q 2 )- l l 2 &q. 
This solution depends on \qo\, and in consequence the re- 
gion of bistability is also |go|-dependent. The bifurcation 
scenario is consistent with two saddle-node (SN) bifurca- 
tions emanating from a (codimcnsion-2) cusp point. 

Case II: \m\/2 < S < \m\; fig. 2(c). The only relevant 
bifurcating lines (in addition to K = \m\) are given by 
K^p in eq. (|25p with a left branch emanating from the 
-ftf-axis and existing up to 7 = 25/|m| — 1, and a right 
branch existing above 7=1. At the left branch of K^ 
the bifurcation from incoherence is supercritical, while it 
is subcritical at the right branch. 
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Fig. 2: Phase diagram for m < with 5 = (a), 0.2|m| (b), 
0.8|m| (c), and 1.2|m| (d). Solid lines are analytical predic- 
tions (tested by numerical simulations). The dashed lines are 
obtained directly from numerical simulations with m = — 1, 
qo = I and N = 2000. The dotted line in panel (d) is a bound 
of the region of bistability given by eq. (|27[) . Small panels show 
magnified regions of the the phase diagrams (a) and (b). 



Case III: S > \m\; fig. 2(d). At 8 = \m\ the locus of 
K^ [eq. (|25p ] reorganizes giving rise to the phase diagram 
for 6 > |to| shown in fig. 2(d). A wedge-like region of un- 
stable incoherence between K^ and K = \m\ exist above 
7 = 2r5/|TO,| — I. This means that in the limit to — > CP 
this region disappears and the phase diagram becomes the 
one found in [T2] in the independent case (m = 0), as ex- 
pected. The fact that the right branch of K^ corresponds 
to a subcritical bifurcation results in a region of bistability. 
This region cannot be analytically determined, though a 
lower bound for its upper border can be calculated finding 
where an unstable solution (with < R < R x ) exists. We 
obtain the line 



Ky 



TO 2 (I - 7) 

2(H-5) 



with 7 > 



2(5 
m| 



1 



(27) 



shown as a dotted line in fig. EJd), which is a good esti- 
mation of the upper border of the bistable region. 

Remarkably, the phase diagrams for negative depen- 
dence differ significantly from those obtained for positive 
dependence (fig. [J). With negative to, synchrony becomes 
more dominant in the phase diagram, in consonance with 



p-4 



The Kuramoto model with distributed shear 



the numerical observations made in [14.. 

Linear stability analysis. — With general distribu- 
tions the residue's theorem cannot be used. We can never- 
theless follow Strogatz and Mirollo 23 and calculate the 
linear stability threshold of the incoherent state. This al- 
lows to know how much our results for the Lorentzian h(q) 
are applicable to other distributions, what is always a con- 
cern when applying the OA theory [13l[21]. Our analysis 
is not completely rigorous but permits to understand the 
results of the numerical simulations. 

In the incoherent state, all Fourier modes (save /o) van- 
ish: fi^ — 0. Equations J7} for the Fourier modes indi- 
cate that at the lowest order only the first Fourier mode 
/ = ±1 is relevant, and it obeys: 

dfx 



(a) 



"-g= - i( W + qK)h 
+ y(l + i<?) 



(28) 

fi(u',q',t)p(u)',q)du]' dq' 



The linear operator in the right hand side has a linear 
spectrum of eigenvalues A. If fi(ui,q,i) = b(ui,q) exp(At) 
is inserted into eq. (f2"5|) and the trivial solution b = is 
discarded, we get: 



2 

K 



1+iq 



I , .. ^ P (oJ,q)dojdq (29) 
l-oo A + i(a; + qK) 

Defining A = A,, + iA^ , one finds that the imaginary part 
of eq. ([2"9"| has always a solution A; = — ujq at the stability 
threshold (A,. — > + ) if the distribution h(q) is centred at 
zero. As an important example, let us mention the case 
of Gaussian PDFs: 

h{q) = j= e~^, g c (u\q) - 



^ 



2ttv 



(hereafter we take cj = as this can always be achieved 
going into a rotating framework). We obtain an equation 
for the stationary (Aj = 0) instability of incoherence: 

K s c \J 2v 2 (K s c +m) 2 + 2a 2 v 2 (K« + m) 2 + a 2 y ' 

For to = a simple analytic solution for K s c can be found 
[T3] ; otherwise K* is the solution of a fourth-order polyno- 
mial. In the next section we show that sometimes complex 
eigenvalues (A-j ^ 0) may also destabilize incoherence, and 
hence using eq. (f30]) we take the risk of missing nonsta- 
tionary instabilities. 

Linear dependence between ui and q. — If there 
exists a purely linear dependence of ui on q, ujj = m(qj — 
qo), general expressions for the stability threshold of inco- 
herence can be obtained if qo = 0. With this latter choice 
the system possesses reflection symmetry (6j , uij ,qj) — > 
{—9j 1 —ujj,—qj) 1 in addition to the rotational symmetry 
6j — > 6j + 4>. Inserting the pdf p(u),q) = h(q)S(u) ~ mq) 
into eq. (|29|) . and taking the limit A r — > + , we obtain: 



2 



7i7t(0) 
\K? + rr 



+ 



1 



(31) 
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Fig. 3: Phase diagrams of model J2} with Gaussian h(q) and 
g c — 5(uj — mq) for m > (a) and m < (b). Solid lines 
(numerically tested) correspond to eq. (|32|l [and K = — m in 
(a)]. Dashed lines are obtained from numerical simulations 
with |m| = 1 and iV = 4000. Small panels show magnified 
regions of the the phase diagrams (a) and (b). 



Solving this equation for K s c gives the boundaries: 



KS = {^m-x if a- > -to, 



(32) 



The linear stability analysis permits to determine at which 
side of the bifurcation the incoherent state is unstable. 
This is an indirect indication that there must exist a 
horizontal bifurcation line at K = —m, exactly like in 
figs. [Ha) and[UJa) for Lorentzian h(q). Moreover eq. ([321) 
agrees with the analytical and numerical results obtained 
in figs.QJa) and[2{a) for positive and negative m, respec- 
tively. Note also that reflection symmetry makes the sta- 
tionary instability at K s c to be a (circle-)pitchfork bifurca- 
tion. We also report next the results obtained with Gaus- 
sian h(q): 

Positive m. The result of our numerical simulations 
with Gaussian h(q) and to > is presented in fig. 3(a), 
and confirms the soundness of eq. ([3"2"|) . In contrast to the 
case of Lorentzian h(q) in fig. HJa), a region of bistability 
between synchronization and incoherence exists at large 
K. This can be understood taking the limit K — > oo, in 
eq. ([2]) and performing a self-consistence analysis a la Ku- 
ramoto, see [13]. A solution branches off from incoherence 
at nh(0) = 1 increasing the value of is, a scenario of sub- 
critical bifurcation coherent with the observed bistability. 
The orientation of this branch is intrinsic to the form of 
h(q) and is independent of the value of to. For distribu- 
tions with a sharp peak, like the triangular or Laplace dis- 
tributions, the bifurcation is supercritical [and the phase 
diagram will be slightly different from that in fig. EJa)]. 
The Lorentzian distribution is marginal and finite- AT ef- 
fects make the bifurcation to be supercritical for to > 
and subcritical for to < 0. 
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Negative m. The numerical results for Gaussian h(q), 
shown in fig. [3Jb) , indicate that eq. ([52]) predicts every- 
where the correct boundaries for stable incoherence, ex- 
cept in a region close to K = \m\ (see bottom panel). 
There incoherence undergoes a Hopf bifurcation at K%, 
a bifurcation line that emanates from a double zero 
eigenvalue (Takens-Bogdanov) point located on K s c at 
vtb — (This stems from a degeneracy at nh(0) = 

— J h'(q)q~ 1 dq.) It is remarkable that the Hopf bifurca- 
tion gives rise to a standing wave (SW) consisting of two 
counter-rotating clusters of locked oscillators. In the stan- 
dard Kuramoto model the SW cannot arise in unimodal 
distributions of w, but it is typical of bimodal distributions 
with well separated peaks 5,19, 20 . The other lines in the 
bottom panel of fig. 131b) are (twin) saddle-node bifurca- 
tions (SN) emanating from a degenerate-pitchfork point, 
and a heteroclinic connection (Het) born at TB. 

Taking <Zo 7^ breaks the reflection symmetry and the 
phase diagram should exhibit structures already found in 
the Kuramoto model with bimodal non-symmetric distri- 
bution [21] or unbalanced interacting subpopulations [25] . 

Conclusions. — Our work is a natural step in the de- 
velopment, initiated by Winfree and Kuramoto, of realistic 
solvable phase models, as simplifications of the mean-field 
complex Ginzburg-Landau equation [BJIH] or in other set- 
ups [1,26,27 . The model analysed in this Letter widens 
the scope of the Kuramoto model by admiting shear diver- 
sity. Shear is a generic feature of oscillators with particular 
relevance in ensembles of limit-cycles close to collision with 
a saddle point (saddle- loop bifurcation) [35]. These sys- 
tems may be good candidates to observe the phenomena 
reported here. 

Considering a broad but still reasonably simple family of 
joint distributions p(u>, q), we have found that the sign and 
magnitude of m, controlling the dependence between the 
natural frequencies and the shears, has a profound impact 
on the phase diagrams. Synchronization is prevalent for 
negative m, whereas incoherence prevails if m is positive 
(or zero |13]L A certainly interesting line of future work 
would be to investigate the effect of other dependences 
between uj and q on the synchronization phase diagrams. 

Finally, this work can also give hints about the valid- 
ity of the OA ansatz in systems with distributed parame- 
ters [25]. Why distributing q is so amenable to analysis? 
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